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Abstract 

One of the most remarkable examples of emergent quasi-particles, is that of the "fractionalization" 
of magnetic dipoles in the low energy configurations of materials known as "spin ice", into free and 
unconfmed magnetic monopoles interacting via Coulomb's 1/r law [Castelnovo et. al, Nature, 451, 
42-45 (2008)]. Recent experiments have shown that a Coulomb gas of magnetic charges really does 
exist at low temperature in these materials and this discovery provides a new perspective on otherwise 
largely inaccessible phenomenology. In this paper, after a review of the different spin ice models, we 
present detailed results describing the diffusive dynamics of monopole particles starting both from 
the dipolar spin ice model and directly from a Coulomb gas within the grand canonical ensemble. 
The diffusive quasi-particle dynamics of real spin ice materials within "quantum tunneling" regime 
is modeled with Metropolis dynamics, with the particles constrained to move along an underlying 
network of oriented paths, which are classical analogues of the Dirac strings connecting pairs of Dirac 
monopoles. 



1 Introduction 

Spin ice materials]!, 2] form part of a series of rare-earth oxide insulator R2M2O7 with space group 
Fdim where R 3+ is a magnetic (Dy 3+ and Ho 3+ ) and M 4+ a non-magnetic (Ti 4+ or Sn 4+ ) ion. The 
cations sit on the vertices of two interpenetrating pyrochlore lattices formed by corner-sharing tetrahedra 
(see figure 1). The well-established members of this family are Dy2Ti207 , H02T12O7 and Ho2Sn207, 
the latter being less studied as it is only available in polycrystalline form. Other potential candidates 
are Dy2Sn207 or Pr2Sn2 07, while related members of the series with rather different properties include 
Tb2Ti207 and Tb2Sn207. The total angular momentum J = L + S is a good quantum number for rare 
earth elements and the electronic ground states are determined by Hund's rules. 

The free ion Dy 3+ (resp. Ho 3+ ) then has a 16-fold (resp. 17-fold) degeneracy that is lifted by the 
surrounding crystal field; the corresponding energy levels have been estimated by neutron time-of-ffight 
spectroscopy [3], revealing an almost pure ground state doublet. 

As the first excited states are separated from the doublet by an energy scale of A s» 2 — 300K (see 
Table 2.), the low temperature behaviour of these materials can be approximated by classical Ising spins 
with a very large magnetic moment (« 10 (J-b), oriented along the line connecting the centers of two 
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Table 1: Ground state electronic configuration of the atom [R={Dy,Ho}], the free ion [R 3+ ], and 
the corresponding quantum numbers S, L and J 
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Table 2: Crystal field level: The ground states | J, J z ) (GS), the Lande factor gj, the estimated magnetic 
moment fi = gj J fj,s and the energy level of the first excited state A. (^gj = 1 + J ( J + 1 )~^(^+i)+ s ( s + 1 ) ^ 



neighbouring tetrahedra and forming a local set of body centered cubic axes: 
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CTj = ±1 is an Ising pseudo-spin with chosen convention that <7j = 1 corresponds to the moment pointing 
out of a down tetrahedron (sec figure 1). 

The physics of spin ice is captured to a good approximation by an effective model with nearest 
neighbour ferromagnetic interactions between the moments which, together with the strong crystal field, 
gives rise to a frustrated ferromagnetic system: 



H = -3 J of f S» • = Jog a i a j- 

(id) (id) 
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Harris & al.[l] showed that the low energy phase space of this system is identical to the Pauling model 
for ice [4] , which shows extensive ground state entropy in which states satisfy the Bernal- Fowler ice rules 
and which has been measured, both in ice ([5]) and spin ice ([6]). The nearest neighbour spin ice (NNSI) 
Hamiltonian also maps onto an Ising antiferromagnet for the pseudo spins [7, 8] (see equation (2)); a 
model first derived by Anderson to describe spinel materials [9] . The equivalence between these three 
nearest neighbour models is illustrated in figure 2. 

The Bernal - Fowler ice rules requiring two magnetic moments to enter and two to leave each tetrahe- 
dron constitute a topological constraint, with the result that the Pauling states make up an ensemble of 
gauge invariant topological sectors with U(l) symmetry. A consequence of such constraints in this and in 
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Figure 1: Pyrochlore lattice: The spins are located on the corner of every tetrahedra and are fixed 
along their local [111] axis represented by the dashed lines. There are two types of tetrahedra that we shall 
call the down tetrahedra (bottom left with the four spins) and the up tetrahedra. Each down tetrahedra 
is connected to four up ones, and vice-versa. The cube represents a unit cell with 8 tetrahedra (four of 
each kind) and 16 spins, and defines the [100] (x), [010] (y) and [001] (z) axes. We introduce the length 
of the unit cell a « 10 A as well as the distance between nearest neighbour r nn = ^ a w 3.5 A and 

between the centres of two connected tetrahedra r c i = ^ a 4.3 A. The smallest close loop encompasses 
6 spins (see green dotted loop). 

related systems is the development of algebraic correlations [10, 11, 12, 13, 14]: consider the magnetisation 
of a tetrahedron a: M(r Q ) = J2iea &i anc ^ define a coarse-grained field M(r) = p. J2aev M(r Q ) averaged 
over a volume V large enough to make M(r) a smoothly varying function. The proof of algebraic correla- 
tions requires a two-stage argument based on the frustrated 2 in - 2 out ground states. On this manifold 
of states, the Gibbs free energy G is purely entropic and can be expressed as a functional of {M(r)}. For 
a tetrahedron respecting the ice-rules, M(r a ) can only take 6 values proportional to (±1,0,0), (0, ±1,0) 
or (0, 0, ±1). Hence computing the entropy density for a volume V is similar to the problem of a random 
polymer made of monomers fixed on a cubic lattice and whose extremities are separated by M(r). A 
large value of |M(r)| requires most of the spins/monomers to be aligned along this direction, allowing 
only a relatively small number of configurations. On the other hand, zero magnetisation fixes a minimum 
number of degrees of freedom. Hence the Gibbs free energy and its Fourier transform are expressed to 
lowest order and up to a constant as 
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Figure 2: Models equivalent to nearest neighbour spin ice: The different mappings between 
Anderson's model for antiferromagnctic spinels, spin ice and water ice: the white spheres are the hydrogen 
atoms and the red one is the oxygen. A spin pointing inside (resp. outside) the tetrahedron corresponds 
to a down (resp. up) spin in the Anderson model and a short covalent bond (resp. long ii-bond) for 
water ice. 

which gives a temperature independent Gaussian probability e~ G l T and (Af Al (— k)M„(k)) = 
However this argument neglects the non-local influence of the ice-rules. If we regard the spins as local 
fluxes of magnetisation then these constraints impose a flux conservation for each tetrahedron (2 spins 
in and 2 spins out)] from a coarse-grained point of view, it is equivalent to a divergence free condition. 

V-M(r)=0 & k-M(k)=0 (5) 

which imposes M(k) to be orthogonal to k and thus 

(M M (-k)M„(k)) = 1 (V - ^) (6) 

where the Greek indices fi, v are Cartesian coordinates labels. After a final Fourier transform back to 
real space, we obtain the desired correlations 

(M M (0)M,(r))ocl ^~ 3 3 ^ where f „ = (7) 

This result deserves a few comments. First of all, we have lost the details of the pyrochlore lattice as 
soon as we used a coarse-grained field and this result can be generalised to models others than spin ice 
[10, 11, 12, 13, 14], the only hypothesis being 

• a Coulomb phase gauge theory imposed by equation (4); 

• a divergence free condition, consequence of the ice-rules constraints (5). 

The magnetisation field M(r) is then analagous to a magnetic field without monopoles which explains the 
specific form of the correlations similar to those generated by a dipole-dipole interaction. One might think 
that the 1/r 3 dependence, in a 3d crystal should lead to a logarithmic divergence in the susceptibility and 
logarithmic peaks in scattering function. This is not, however the case as the dipolar angular dependance 
of the correlation function results in a cancellation of the divergent part when integrated over all space. 
The correlations can be observed experimentally however through the development of narrowing pinch 
points [15] in the diffuse scattering profile [16, 17, 18]. 

It is clear that the NNSI model already contains extremely rich phenomenology related to the con- 
strained Coulomb phase with U(X) gauge symmetry, which is the basis for monopole physics. Most 
notably, single spin flip excitations create a pair of topological defects; pairs of tetrahedra with 3 in - 
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1 out and 3 out - 1 in . The defects are expected to have an effective entropic interaction [14], but 
once created can move at zero energy cost. They therefore posses some but not all the characteristics of 
neutral magnetic monopole pairs. In order to see the emergence of true magnetic quasi-particles one has 
to go beyond the NNSI through the inclusion of dipolar interactions, as is considered in the next section. 



2 The dipolar spin ice model 

Dipolar interactions are of particular importance [19, 20] for the rare earth ions Ho 3+ and Dy 3+ at the 
heart of spin ice physics: as the exchange coupling between the ions is due to 4/ electrons, buried behind 
5s and 5p layers 1 , it is very weak (~ IK) in comparison with other ferromagnets (for example J <~ 10 2 
K for metallic iron). The magnetic moments in spin ice being very large, the dipolar interactions, 
that are negligible on short length scale for iron, are on the same energy scale as the exchange here. 
A quantitative description has therefore to take both interactions into account while to an excellent 
approximation moments can be described by Ising degrees of freedom orientated along the body centers 
of the tetrahedtra. This separation of scales leads to the dipolar spin ice (DSI) Hamiltonian [20] 
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where J, D and r nn w 3.5 A are respectively an antiferromagnetic exchange coupling, the dipole-dipole 
coupling and the nearest neighbour distance between rare earth ions, which has been shown to provide 
a comprehensive quantitative description of spin ice materials. Here, J is a free parameter, while D is 
defined through the dipolar interaction 
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for Dy 2 Ti 2 07. A numerical estimate of J has been obtained by comparison between Monte Carlo simula- 
tions and experimental data for the specific heat in Dy 2 Ti 2 07 [20] or of neutron scattering on Ho 2 Ti 2 07 
[21] and Ho 2 Sn 2 C>7 [22]. Uncertainties in the estimated values for the Stanates remain high; around 
50% for Ho 2 Sn 2 07 because of the unavailability of single crystals. The dipolar contribution turns out to 
be indispensable for understanding certain magnetic field induced transitions [23, 24, 25, 26] as well as 
the quantitative feature of magnetic monopole excitations [27]. However, a remarkable property of this 
Hamiltonian is that the extensive degeneracy of an effective nearest neighbour model is almost main- 
tained, as the nearest neighbour contribution to the dipolar interaction is the same for all Pauling states 
and the long ranged part of the interaction is almost perfectly screened[28, 29]. The dipolar interaction 
between neighbouring spins 1 and 2 defined in figure 2 is given by: 
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By symmetry, the result is the same for all pairs of spins in a tetrahedron. Hence the NN term of 
the dipolar interaction is equivalent to an exchange term and the nearest neighbour spin ice (NNSI) 
Hamiltonian is recovered as follows 
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Table 3: Nearest neighbour interactions: Exchange, dipolar and effective anti-ferromagnetic cou- 
plings in Kelvin. 

Values for the spin ice materials is summarised in table 3. below. 

The origin of this approximate screening, the so called projective equivalence, [29] can be seen phys- 
ically by considering a further abstraction of the DSI to the dumbbell model [30], which highlights the 
remarkable symmetry properties of dipoles on the pyrochlore lattice. In this model the point dipoles of 
the DSI are extended into magetically charged dumbbells, as shown in figure (3), whose ends sit on the 
centers of the tetrahedra, forming a diamond lattice of magnetically neutral sites for all Pauling states. 
Projective equivalence is clearly exact for this model and the error incurred in DSI, which turn out to be 
quadrupolar corrections only [29] is due to the difference between point and extended dipoles. 




Figure 3: Dumbbell model: Each magnetic dipole is seen as two charges sitting on the vertices of the 
diamond lattice (dashed lines). 2 in - 2 out is a vacuum whereas 3 in - 1 out is a positive charge. 



3 Magnetic monopoles and classical "Dirac strings" 

Within the dumbbell model, the ensemble of Pauling states provide a quasi-particle vacuum for neutral 
pairs of topological defects. Once created, the 3 in - 1 out and 3 out - 1 in defects break magnetic 
charge neutrality on each duel lattice, so that, when separated, the defects should interact via Coulombs 
law. Hence, one can see by construction that fractionalization of the dipolar dumbbells into deconfined 
magnetic Coulomb charges occurs when the ice rules are broken: in this sense these quasi-particles are 
magnetic monopoles. To show that, thanks to projective equivalence, the same is true to an excellent 

lr The 6s layer is empty because Ho 3 " 1 " and Dy 3 + lost three electrons. 
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approximation in the DSI, we have followed [30] and simulated a DSI system. The long range interactions 
have been treated using the Ewald method [31, 32, 23], where the slowly converging sum of all dipolar 
interactions between spins is replaced by two strongly convergent series. Periodic boundary conditions 
ensures correct convergence. Starting with a random 2 in - 2 out configuration, we can create and force 
the diffusion of a single pair of monopoles by flipping spins and then compute the energy of the system 
for each configuration. After averaging over a large number of initial microstates and paths of diffusion, 
we obtain the potential of interaction between two monopoles of opposite charges shown in figure 5. For 
a distance r — xr^ where r nn is the distance between two vertices on the diamond lattice (see figure 1) , 
the Coulomb potential can be written [30] 

t t" / \ Mo Q 2 V min 2[i 

V(r) = — = , where Q = — , (12) 

At: r x r<i 

and where 

v - Mo 4 M 2 1 _ _8 

with, from equation (9), D w 1.41 K. V m - ln is the energy gained, with respect to the vacuum, by the 
creation of a nearest neighbour pair of monopoles. This potential V for the dumbbell model is compared 
with our numerical results for DSI in figure 5. As we have periodic boundaries, the pair of monopoles 
under consideration will also feel the presence of its image charges separated by at least one system 
size. While this is convenient for simulating a real system (as is done later in this paper), it leads to 
a finite size bias to the Coulomb energy here, where we consider a single pair of monopoles. This bias 
remains negligible for small x but can be observed at larger x and scales away with system size. This is 
illustrated in 5 for two system sizes L = 4 (green crosses) and L = 8 (red dots). For the larger system 
size the 1/x behaviour is clearly observed and the data respects the energy scale fixed by V m i n to a 
good approximation. The data plotted here are averaged over at least 100 simulations, so that the small 
differences between our results and Coulomb's law reflect the limits of accuracy of the dumbbell model, 
which arc of order 1/r 5 . The only fitting parameter in the figure is the reference energy scale for the 
DSI, which translates the numerical data along the axis. Similar results are obtained in reference [30]. 
However while they computed the energy of separation between monopoles for a single configuration, we 
have performed a statistical average, which produces some statistical noise, but confirms that this is a 
general property of the ensemble of constrained Pauling states. 
Table 4. summarises the main parallels between the dumbbell and DSI models. 
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Table 4: Mapping from the dipolar spin ice to the dumbbell model 
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Figure 4: Emergence of monopoles a) The magnetic ions (Ho 3+ or Dy 3+ ) lie on the sites of the 
pyrochlore lattice and are constrained to the bonds of the diamond lattice (dashed lines) . Local topological 
excitations 3 in - 1 out or 3 out - 1 in correspond to magnetic monopoles with positive (blue sphere) or 
negative (red sphere) charges respectively, b) The diamond lattice provides the skeleton for the network 
of Dirac strings with the position of the monopoles restricted to the vertices. The orientation of the Dirac 
strings shows the direction of the local field lines in H. 

The theoretical basis for magnetic monopoles in the DSI model is thus clearly established. We 
particularly stress that the scale of the Coulomb interaction is set by magnetic constants: the permeability 
of free space, /j,q, and the nearest neighbour distance on the diamond lattice. Hence the quasi-particlcs 
behave quantitatively as magnetically charged particles. However, two fundamental differences still exists 
between these particles and Dirac monopoles: firstly they correspond to divergences in the magnetic 
intensity H, or magnetic moment M, rather than in the magnetic induction, so that V B = V-(H + M) = 
and they do not require a modification of Maxwell's equations. Secondly, their magnetic charge is not 
quantified. The magnetic charge on a Dirac monopole is quantified through a quantum conjugation 
relation with the fundamental electronic charge [33, 34]. This is clearly not the case in spin ice and in 
fact the charge can be continuously varied by modifying the nearest neighbour distance between spins by 
applying external pressure [30]. However, on all length scales above the atomic scale, a 3 in - 1 out defect 
appears to be a local sink in M and therefore a source of field lines in H. Furthermore we can assign 
to them a positive or negative charge when immersed in a magnetic field [35] that is conserved through 
the creation and annihilation of neutral monopole pairs. As the Coulomb potential is not confining 2 in 
3d, these quasi-particles provide an example of fractionalisation in high dimensions: the point dipoles of 
the DSI appear to be separated into free poles of opposite sign that are free to wander independently 
through the system. They are therefor classical analogues of Dirac monopoles and spin ice provides the 
first 3d experimental realisation of deconfined magnetic monopoles. 

In the theory of Dirac, monopoles of a neutral pair of charges are connected by a tensionless Dirac 
string of overturned dipoles, taken in the the continuum limit. The string, which corresponds to nodes 
in the wave function for the pair of particles, should be unobservable. Classical analogues of Dirac 
strings exist in DSI and in the dumbbell model [30], which correspond to the magnetic moments in the 

2 having a Coulomb gas on a lattice prevents the divergence for r — > + . The energy of deconfinement is finite and equals 
to |V m i n | on average. 



8 



X 

> 




Figure 5: Coulomb potential: Effective interaction (in Kelvin) between a pair of monopoles in the 
dipolar spin ice model as a function of the distance between them in units of nearest neighbour distance 
between monopoles. The system sizes are 4x4x4 unit cells, i.e. 1024 spins (green crosses) and 
8x8x8 unit cells, i.e. 8192 spins (red dots). The solid line is the Coulomb potential expected from 
the dumbbell model V(x) — V m i n /x. As we are on the diamond lattice, the smallest possible distance 
between monopoles is x = 1. The data are the average of over more than 100 samples. 



constrained 2 in - 2 out states. That is, by creating and separating a pair of monopoles one flips a series 
of moments, defining one possible Dirac string connecting them. Starting from a disordered 2 in - 2 out 
state, the path taken to create the monopoles and place them at positions A and B is not unique and 
the arbitrariness provides a classical analogue of the quantum uncertainty of Dirac monopole trajectoires 
in space. Further to this analogy, within the dumbbell model, as long as the dumbbells are taken to 
have point magnetic charges at each end (the dumbbells can also be thought of as infinitesimally thin 
magnetic needles, touching at the centres of the tetrahedra), the magnetic flux associate with each dipole 
is confined to infinitely narrow elements joining the charges on each dumbbell, so that the network of 
classical Dirac strings is also invisible. This of course is a mathematical abstraction; in the DSI, we have 
point dipoles on the vertices of the tetrahedra and because of this difference (because of corrections to 
projective equivalence in fact), the string network is observable. Indeed, recent experiments have reported 
the observation of "Dirac strings" [17, 36]. As far as dynamics are concerned, the Dirac string network 
plays a crucial role: moving a monopole involves a spin flip, so that the quasi-particles leave a wake of 
re-arranged links in the network of Dirac strings as they move, with 3 in - 1 out and 3 out - 1 in defects 
moving in opposite directions over the network (see figure 4). This directionality makes the dynamics 
highly constrained compared with the diffusing electric charge of an equivalent electronic Coulomb gas. 
At equilibrium and in zero field this constraint simply renormalizes the diffusion constant but in a 
rapid quench the constraints lead to glassy behaviour [37] reminiscent of other kinetically constrained 
models[38]. We will return to the dynamics in finite field in the conclusion. 

Given the accessibility of these magnetic quasi-particles, the development of an experimental signature 
is of vital importance and interest and has recently given rise to considerable experimental activity[17, 36, 
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18, 39]. In particular, recent experiments have measured the charge carried on free deconfined particles 
via the magnetic Wien effect [40, 41]. Here, we will concentrate on magnetic relaxation, available from 
measurements of magnetic response in weak ac magnetic helds [42, 43, 44, 45]. In the following section 
we will show how these data can be used as a direct signature for diffusive monopole dynamics. 

4 Dynamics in spin ice 

In figure (6) we show the evolution of the magnetic relaxation time scale for Dy 2 Ti 2 07 measured from bulk 
a.c. susceptibility measurement [45]. The data show three distinct features: a high temperature regime 
above 12 K, which, has been associated with thermal excitations of higher crystal field levels [46, 44, 47], 
an intermediate temperature plateau regime between 3 and 10K and a regime with rapidly increasing 
time scale below 3 K. Qualitatively similar behaviour has been observed in Holmium Titanate[48]. We 
will use this data to illustrate that the low and intermediate temperature dynamics of spin ice can be 
interpreted in terms of the stochastic dynamics of monopole quasi-particles. Given the separation of 
energy scales in spin ice materials we follow the assumption that the dynamics in the low temperature 
and the plateau regions are controlled by quantum tunneling events [44] taking the magnetic moments 
directly from one Ising position to another through the existence of non-zero off-diagonal components of 
the dipolar interaction from neighbouring spins, giving a small mixing of the magnetic states. Within this 
picture an excitation requires thermal assistance to make an Ising spin flip, but the time scale for the flip 
should be given by rate constant t , which is essentially temperature independent. The high temperature 
limit of this model is then an Ising paramagnetic system evolving with temperature independent rate 
constant, giving a plateau in the evolution of the relaxation time scale. The deviation from the plateaux 
above 10 K is due to thermal mixing of the crystal field levels [46, 44, 49, 50, 51, 52] and will not be 
discussed further here. 

From figure (6) the tunneling rate, r for Dy 2 Ti 2 07, appears to be around 10~ 3 s , while for Ho 2 Ti 2 C>7, 
it is on the scale of r Q w 10~ 8 s [48]. These extremely slow, yet exponentially varying microscopic time 
scales from one compound to another, seem consistent with a quantum tunneling process exponentially 
dependent on the energy levels of the rare-earth ions CEF. Unfortunately, the complexity of these CEF 
prevents quantitative estimates for the experimental inversion rate and the time scale t q remains a 
fitting parameter in our theory. The quantum tunneling on these microscopic time scales should lead to 
stochastic dynamics on all length scales and if a single micropscopic process dominates the dynamics, one 
should be able to model the relaxation using an Arhenius law. For a more complex system, for example 
one with long range interactions, one should be able to use the Metropolis dynamics of a Monte Carlo 
simulation to model the evolution of the system. These techniques and the validity of the model, with 
temperature independent r are discussed in detail in the next section. 

4.1 Arrhenius law and deconfined quasi-particles 

The first quantitative analysis of spin ice dynamics [42] relates to Ho 2 Sn 2 07 and Ho 2 Ti 2 07. Here the 
relaxation time, t(T) was fitted with an Arrhenius law in the freezing regime (T w 1 K) with characteristic 
energy barrier, Ej — 19.6 K and 27.5 K respectively. These numbers already suggest that a simple 
Arrhenius law cannot give a complete description of the dynamics, as they do not correspond to any of 
the energy scales appearing in the DSI: 

• ~ 200 K between the ground state doublet and the first excited crystal field levels; 

• ~ 4J cff = 7.2 K for Ho 2 Ti 2 7 and 4.4 K for Dy 2 Ti 2 7 , due to single spin flips within the NNSI 
model; 

• ~ 3 K for the limit of infinite separation between monopoles (see figure 5). 
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Figure 6: Experimental relaxation time r for Dy2Ti207 vs temperature T from Snyder & al. [45], 
courtesy of Peter Schiffer. The starting point of this work has been to shed new light on this result. 
Above 12 K the spins are not exactly Ising anymore. The quasi-plateau between 2 and 12 K is divided in 
two parts: the "high temperature" of the spin ice region above 4 K where one needs to consider double 
defects 4 in or 4 out and a large monopole density, and the "low temperature" regime below 6 K where 
our monopole model should hold. Below ~2K, the system becomes frozen. 



Quantitative studies of Dy 2 Ti 2 07 rule out such a simple thermally activated process [43, 45], as shown 
in figure 7. Here, the data of Snyder & al. arc fitted in the low temperature regime with an energy 
barrier of Ef — 6J c s ~ 6.6 K, with J c s = 1.11 K, the value estimated for Dy 2 Ti 2 07 [20]. The Arrhenius 
scaling is quantitatively good below 2 K, but fails completely to reproduce the quasi-plateau region. As 
we have argued above, if the tunneling picture is correct, the plateau region, through its simplicity, offers 
an anchor point for any fit and any procedure that does not fit here cannot contain the physics of the 
problem. 

Rather than fitting the low temperature freezing region one should therefore start by fitting the 
data through the plateau region. By doing this we see the first evidence for deconfined quasi-particle 
excitations. Following the idea that the quasi-plateau is in fact the tail of an exponential, thermally 
activated process, we have plotted in figure 7 different Arrhenius laws To exp(_E p /T) where the energy 
barrier is varied. The time scale To is fixed by scaling to the experimental time at 4 K; a temperature 
at which we expect the density of double defects (4 in and 4 out) to be negligible, leading to a direct 
relationship between spin flips and monopole pair creation. From the figure one can see the main result 
of this section: the experimental data are poorly reproduced by an energy barrier E p = 4 J e g, the cost 
of a single spin flip within the NNSI model, but are quantitatively well fitted with E p = 2 J e g over the 
region between 2.5 K and 5 K (see the lower panel of figure 7). This means that the lowest excitation 
responsible for the dynamics of spin ice materials is not the energy cost of a single spin flip, but only 
half of that i.e. the energy cost of a single topological defect. Within this effective temperature range, 
the NNSI model seems sufficient to capture the physics, as suggested by polarized neutron scattering 
[53] . Once a pair of defects is created, they can propagate freely and separately in the system, leading to 
magnetic relaxation from independent (deconfined) objects. 
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Figure 7: Deconfined monopoles: The experimental data of Snyder & al. [45] (x and dashed line) are 
compared with different Arrhenius functions. Top: The quasi-plateau region is in quantitative agreement 
with a thermally activated process with energy barrier E p — 2 J c g (red line) whereas the spin freezing 
is well reproduced with E p rj 6 J e g (green line) , but no unique function can fit the whole temperature 
window. Bottom: The characteristic excitation is the creation of a unique defect (E p = 2 J e g, red line) 
rather than a single spin flip (E p = 4 J e g, grey line). The fit is improved at higher temperature if we 
include the energy scales due to double defects in a modified Arrhenius law (blue line). 
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As shown in figure 7, some corrections to Arrhcnious behaviour develop above 5 K, which we show in 
the next section, is due to the excitation of 4 out and 4 in double defects, by constructing a variant of 
the Arrhenius law, including all relevant energy scales for single spin flips within the context of the NNSI 
model. 



4.2 Multi-energy Arrhenius law 

We use an independent tetrahedron picture, assuming that the probability of finding a tetrahedron in a 
certain state is simply given by the Boltzmann weights of the two tetrahedra involved in any spin flip. 
We label the different configurations 2 in - 2 out , 3 in - 1 out (or 3 out - 1 in) and 4 in (or 4 out) 2, 
3 and 4 respectively. In table IV. 2 we show a list of all energy changes resulting from a single spin flip 
shared between two tetrahedra, from which we can calculate the mean energy barrier AE. 
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Table 5: Energy scales resulting from a single spin flip: Let us explain this table using the second 
row as an example. It means that if we flip a spin that initially belongs to a tetrahedron 2 in - 2 out and 
a tetrahedron 3 in - 1 out , then there is a probability of 3/4 simply to inverse the position of the defect 
(3 1 2) at no energy cost, and a probability of 1/4 to create an additional pair of defects (3|4) that cost an 
energy SE = +8 J c ff. The probabilities are a statistical average, whether the spin under consideration is 
one of the three in spins of the 3 in - 1 out tetrahedron (5E = 0), or the fourth out spin (SE = 8 J e g). 

There are two points worth noting: the configurations (2|3) and (3|2) are different, which gives a 
factor of 2 for their Boltzmann weight. More importantly, a single spin flip that actually creates a pair 
of defects will effectively only cost one half of SE, because of the fractionalisation of the excitations, but 
the process (3 1 3) — >(2 14) really costs SE = +4 J c ff as it simply moves a charge. We find AE: 

2 (6 e^f + 4 (2 \ 6 8) + 6 (± 8») + 4 ( jL g) 

^ ' la „2B,l.. tt \ 2 , (o 1 a „1BJ.« cA i I 1 q2\ i i 6 o2\ cff 



T)'+(2l6e^8) + (^82) + (&82) 
6 e 4 / 3 J ° !! + 8 e 2f3J ^ + 10 



- u 9e 4/3J ott + 6e 2/37 eff + 7 J * s T^o 2 Joff ' (14) 

from which we recover the creation energy of a single defect in the low temperature limit. The modified 
Arrhenius law To exp(AE(T)/T) is plotted in figure 7 (blue curve). It is almost identical to the thermally 
activated process with E p — 2 J e g (red line) below 4 K, but differs at higher temperature and follows 
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the experimental data more closely up to w 8 K. Quantitative differences develop above this limit as one 
expects: above this temperature the thermally activated processes of the high temperature regime should 
make their presence felt. Using the function tq exp(AE/T) clearly extends the range over which this phe- 
nomcnological model can describe the data and we believe that it convincingly confirms the validity of the 
quantum tunneling model to describe the dynamics up to 8 — 10K. We believe that this phenomenologi- 
cal approach provides a very strong argument in favour of magnetic relaxation assisted by the thermally 
activated creation of deconfined quasi-particles in spin ice and hence in favour of fractionalisation of 
dipoles into free effective monopoles ! 

However as shown in 7, any Arrhenius function ultimately fails below 2-3 K, underestimating the time 
scale in the low temperature region. Along the quasi-plateau the density of monopoles is high enough for 
Debye screening to occur with the result that effective models based on the NNSI provide a convincing 
description. As the temperature is lowered, their density is reduced and a quantitative description can 
only be achieved by taking explicitly into account the long range Coulomb interactions. In this case, 
simple Arrhenius behaviour is excluded and we must revert to numerical simulation. 

5 Spin freezing 

We have simulated the dynamics of both the DSI model and a Coulomb gas of magnetic monopoles 
in the grand canonical ensemble occupying the sites of the diamond lattice, using a Metropolis Monte 
Carlo algorithm. The Coulomb particles are excitations out of the quasi-particle vacuum provided by the 
dumbbell model and their dynamics is constrained by the underlying network of Dirac strings. In both 
cases we treat the long range interactions, cither dipolar or Coulombic using the Ewald method [31, 32]. 

5.1 Dipolar spin ice 

Simulation of the DSI model is relatively straightforward. We implement both nearest neighbour and 
dipolar interactions, using the numerical values of J and D for Dy2Ti2C>7 given in table 1.3. We extract 
the characteristic relaxation time t, from the auto-correlation function 



where N is the total number of spins (N = 16 x 8 3 = 8192 here) and Sj(i) is the value of the unit vector 
representing the spin at the Monte Carlo time t. One Monte Carlo time step, for a system of size N 
is defined as a single sweep through N moments chosen at random. For the initial conditions we take 
a random configuration, which we let evolve at high temperature T = 10 K until equilibrium, defined 
here as the time beyond which C(t) decays below 0.01. This defines t = and C(t) is then computed 
and stored for a given temperature. The Monte Carlo time is re-set to zero when C{t) decays beyond 
0.001 and the temperature T is lowered by ST. The process is repeated until we reach T = 0.7 K when 
numerical equilibration becomes difficult. The resulting auto-correlation is averaged over 200 samples in 
order to give a smoothly decaying function down to C(t) < 0.03. The outcome is plotted in figure 8 for 
T =1.5 and 3 K. 

C(t) relaxes almost exponentially, but with a small difference between short and long time scales. 
Different time scales, Ti = \^, can be extracted by fitting the correlation function C(t) to an exponential 
function, cither over the interval C = 1 — 0.8 ( n) or over the interval C = 0.3 — 0.03 (r2)(see dashed 
lines in figure 8). The characteristic relaxation time r, is taken to be the average of t\ and T2- The 
Monte Carlo time is again scaled to that measured by Snyder & at by equating the scales at 4 K, 
t(T = 4 K) = 2.99 ms, for which the density of double defects is negligible (< 1%), giving a quantitative 
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Figure 8: Time correlation function C(t) for DSI for T = 1.5 (green) and 3 K (red) on a semi-log 
scale. The x-axis is the Monte Carlo time t. The dashed lines are exponential fits for short and long 
relaxation times. 

equivalence between our Metropolis dynamics and the relaxation dynamics of the experiment, with one 
Monte Carlo step = 2.5 ± 0.5 ms. 

Preliminary results (not shown) were for the NNSI model using the J c g from Table 3. The data agree 
quantitatively with the analytic calculation, To exp(A_E(T)/T), (see figure (7)) over the temperature 
range — 10 K showing that our rather simple approach accurately describes the stochastic dynamics 
of the NNSI model. Any deviation from this behaviour is therefore explicitly due to the long range 
interactions of the DSI. Our results for r vs T for the DSI are compared with the data of Snyder & al. 
[45] in figure (9). We find remarkably good agreement between the numerical results and experiment over 
the temperature range between 1 K and 10 K. In particular there is a substantial improvement over the 
Arrhenius law of figure (7) below 3 K as we move into the low temperature regime with rapidly increasing 
time scale. Although when one looks more closely there are still some systematic differences between the 
data from the real and the numerical experiment, we believe that the general agreement shown in figure 
(9) represents a considerable success for both the DSI and stochastic quantum tunneling hypothesis, in 
which the tunneling rate is temperature independent. Further, the tunneling time scale is approximately 
equal to the Monte carlo time scale fixed at 4 K, giving To ~ 2.5 miliseconds. 

We start the investigation of exactly how the Coulomb interactions modify the spin ice dynamics by 
calculating the evolution of topological defects as a function of temperature in both the NNSI and DSI 
models. The evolution is shown in figure 10. In both cases the density of 3 in - 1 out and 3 out - 1 in 
defects shoots rapidly to zero below 1 — 2 K, which corresponds to the low temperature regime, while 
the double defect concentration goes to zero below the 5 — 6 K range. This confirms the picture that 
the plateau region is due to the proliferation of a high density gas, while the low temperature region 
with increasing time scale is driven by an exponential fall in the monopole concentration. It is therefore 
clear from this picture that monopole physics should explain the agreement between real and numerical 
experiments shown in figure (9). A second point to notice is that the monopole density falls to zero at 
a noticeable higher temperature for the DSI than for the NNSI model which must surely be related to 
the brutal slowing down of the dynamics compared with our analytical Arrhenius law calculation. To 
investigate this slowing down further we move to the dumbbell model and a description of the problem 
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Figure 9: Relaxation time scale from DSI The experimental data (x) are from Snyder & al. [45] 
(•). The system size is N = 16000 spins. The details of the simulations are given in the text. The 
bars do not represent the error bars on the value of r; rather they show the extreme values T\ and T2 at 
each temperature. In order to see the improvement due to dipolar interactions, the modified Arrhenius 
formula 14 is also plotted in blue. 

uniquely in terms of monopoles. 

5.2 Monopoles and strings 

Within the monopole picture spin flips out of and into the constrained manifold correspond to cre- 
ation/annihilation of pairs, while all other flips correspond to monopole hopping between nearest neigh- 
bour sites of the diamond lattice. Hence in a Monte Carlo step one considers two nearest neighbour sites 
at time t. The possible outcomes at time t + 1 are: 

• If there are no quasi-particles on these two sites, then we consider the creation of a pair of opposite 



• If there is only one monopole, it can either move to the other site or stay where it is; 

• If there are two opposite charges, they can annihilate; 

• If there are two charges of the same sign, no movement is possible; as we have disallowed monopole 
defects carrying double charges. 

5.2.1 Chemical potential [i 

Moving into the monopole picture means we change the Gibbs ensemble from which we select microstates. 
In the magnetic picture the independent variables are T, magnetic field B (chosen to be zero throughout 



charges. 
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Figure 10: Density of monopoles for the dipolar spin ice model (•) and their expected Boltzmann 
weight for the nearest neighbour spin ice model (solid line) with the numerical values of Dy2Ti20y. In 
red (resp. blue) are plotted the concentration of 3 in - 1 out and 3 out - 1 in defects (resp. 4 in and 4 
out), carrying a magnetic charge ±1 (resp. ±2). The system size is N = 16000 spins. 

this paper) and the number of spins N, which is also fixed to the volume of the sample. In the monopole 
picture, the number of quasi particles varies, so the independent variables are T, B — 0, the volume 
V and \x the chemical potential for monopole pair creation. The relevant thermodynamic potential is 
therefore the Grand Potential $ = U — ST — Nfi, for which we have to define \i. We are, of course free 
to choose any value of /j, and arbitrary values will take us away from spin ice physics, but as the present 
work is motivated by the experimental results of Snyder et. al. we need to choose the value of /i that 
best corresponds to Dy2Ti207. In a first series of simulations we have estimated /i = /ii numerically by 
calculating the difference between the Coulomb energy gained by creating a single pair of neighbouring 
magnetic monopoles of opposite sign AZ7 mcmo < in a vacuum and that required to produce a unique pair 
of topological defects out of the 2 in - 2 out manifold in the DSI, AJ7 dof > 0, giving a configurationally 
averaged estimate = Af7 mono — AUdd = —8.92 K. Note that while AU mono is simply V m i n from 
equation (13), AUdef has to be computed by simulation, as projective equivalence is not exact for the 
DSI. In a second series of simulations, fJ-s(T) was taken as the value required to reproduce the same 
defect concentration as in a simulation of dipolar spin ice at temperature T, as shown in figure 10. We 
found that [ii — [i\ within numerical error at T = 0.7K and varied by «3% only, over the temperature 
range 1 — 4 K, showing that our procedures are consistent and underlining the fact that the chemical 
potential used is not a free parameter. Rather, it is taken from detailed comparison with the DSI that 
best describes Dy2Ti2C>7, with parameters from Table 3. Note that in references [27, 54] we have used 
an unconventional definition of fi, with a sign change compared to the above value. Here, we make 
a definition compatible with conventional statistical mechanics such that, in a Monte Carlo move we 
compute the change in SU = 5U mono — ji5N . As the value given, [i w —8.9 K, is for a pair of particles, 
SN = ±1 or zero. 
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It is interesting to dwell further on /j, and its composite parts at finite temperature, as their evolution 
with T turns out to be at the origin of the rapid increase in the time scale at low temperature. We 
can make a third estimate for the chemical potential; n(T) — AU mono (T) — AUdd(T) using the same 
procedure as for fix, but at finite temperature, where a pair of monopoles is created in a microstatc 
containing a finite concentration of monopoles corresponding to finite temperature T. To do this, we 
simulated the DSI model, with TV = 8192 spins in parallel with a 'slave' Coulomb gas, i.e. for each single 
spin flip the corresponding monopole configuration on the diamond lattice is updated accordingly, making 
it possible to compare the energy of the spin configuration with its mirror in the dumbbell model at each 
Monte Carlo step. As we are interested in the creation of monopoles, we have artificially suppressed 
double defects in the DSI for all temperatures. Hence after equilibration at a given temperature T and 
for each accepted creation of a pair of defects, we compute the energy difference between the Monte Carlo 
move in the DSI model AE/dcf and the evolution in the slave Coulomb gas A[/ mono , which gives a value 
of /x for this specific move. By averaging these values over the Monte Carlo time and for different initial 
configurations, one finds the temperature evolution of these energies as well as their standard deviation, 
plotted in figure 11. 
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Figure 11: Creation of a pair of quasi-particles: Energy required in the dipolar spin ice model 
(AUdef , red •) or gained for a Coulomb gas (AU mon o, blue ■) and the resulting chemical potential /i(T) 
(A) as a function of temperature. The vertical bars are not error bars but the standard deviation of these 
quantities. /z(T) tends to the limit \i\ = —8.92 K (lower dashed line) at very low temperature. The other 
dashed line is the analytical prediction of the energy gained by creating a pair of monopoles in a vacuum 
V m i n — —3.07 K (see equation (12)). The temperature scale below 3 K is enlarged for a better display. 



Much of the physics of spin ice at low temperature is included in this figure: the standard deviation 
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of both AUdef(T) (red) and AC/ mono (T) (blue) are large above 1 K, up to ±40% of the mean values. 
However, these quantities are highly correlated and their difference, giving rise to /u(T) (black) is very 
precisely defined with a standard deviation of ±4% at most. Furthermore, the extracted /z(T) follows 
closely the values for /Z2, which was estimated more phenomenologically. This agreement, despite the 
large fluctuations, illustrates the success of the monopole picture and of the dumbbell model in describing 
the physics of spin ice: the energy of separation of two monopoles is given by the magnetic Coulomb 
potential, including many body effects for a finite particle concentration and the energy of creation (and 
annihilation, not shown here) is consistent with a Coulomb gas in the grand canonical ensemble. On 
average, large energy cost AC/d e f (top of the red bars) correspond to small energy gain AC/ mono (top of 
the blue bars) and vice-versa. 

The data in figure (11) show two regimes of behaviour: from 10 K down to 2 K, corresponding to the 
plateau region for relaxation times and below 2 K where the relaxation time rapidly increases. In the 
higher temperature regime the mean values, (AC/d c f) an d (A[/ mon o) are roughly constant and (A£7<jef) = 
(A£7 m0 no) — (M^ 1 )) i s cl° se to the energy cost of a single spin flip in the NNSI model E = 4 J c s = 4.44 
K. Thus, apart from the fact that there is a quasi-continuous range of energy values to create monopoles 
instead of a unique one and that the fluctuations in the energy scales are surprisingly large, the results 
are consistent with the NNSI model, which explains why the NNSI provides a reasonable description of 
spin ice down to 2 K. This equivalence between models with long and short range interactions is a direct 
consequence of Debye screening [36, 14]. That is, throughout this regime, (AC/ mono ) (blue curve) is in 
excess of the energy of deconfincment of a pair of monopoles in a vacuum, V m % n = —3.07 K (lower dashed 
line, sec equation (12)) and from the fluctuations one can see that pair creations providing less energy are 
rare. Hence, a newly created pair is more stable than in a vacuum with the additional energy coming from 
the build up of correlations within the Coulomb gas, favouring a neighbouring charge cloud of opposite 
charge around each particle. As the temperature is reduced, screening becomes more marked so that low 
energy events become increasingly rare, leading a maximum in the mean Coulomb energy gain on pair 
creation, at around 2.5 K. However, below this temperature the monopole concentration drops rapidly 
(see figure 10) so that the Debye screening is starved and there is a marked evolution of the energies in 
figure (11). This appears to occur as the smallest values of AC/ mon o disappear brutally and (AU mono ) 
increases rapidly towards V m in- This evolution is accompanied by a surprising increase in (AC/def) (red 
line), ensuring that fi(T) remains temperature independent. The increase in the mean energy required to 
create a pair of topological defects as temperature decreases is clearly due to the absence of many body 
screening effects from the mobile quasi-particles. The brutal non-Arrhenius spin freezing as one leaves 
the plateau region at low temperature is therefore an avalanche effect: fewer monopoles hinders the 
Debye screening, which reduces the number of creations of quasi-particlcs, which leads to less monopoles, 
... If this slowing down is sharp but continuous, it is because the system remains thermally activated. 
That is, as long as the material remains "hot" enough to allow the creation of a pair of monopoles out of 
the vacuum, i.e. a single spin flip of energy AC/def = Mi + Knm ~ 5.8 K, the system will remain dynamic 
through the creation and diffusion of monopoles. Our simulations hit this limit for T w 0.6 K and one 
might expect equilibration in experiments to become dramatically difficult below this temperature. This 
is exactly what Snyder & al. observed at 0.65 K for Dy2Ti20y in FC-ZFC measurements [45], although 
magnetocaloric measurements [55] and recent muon spin resonance experiments [41] appear to measure 
the dynamical properties of spin ice materials at considerably lower temperature where the intrinsic 
monopole concentration, as predicted by the DSI or the dumbbell model would be considerably lower. 

5.2.2 Dirac strings 

The monopoles hop between nearest neighbour sites via the Metropolis Monte Carlo algorithm, giving 
diffusive dynamics, but with the local constraint discussed in the introduction: in the spin model a 3 in - 
1 out topological defect can move at low energy cost by flipping one of the three in spins, the direction 
of the out spin being barred by an energy barrier of 8 J e g that would result in a double defect and which 
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is set to infinity in the Coulomb gas. An isolated monopole can therefore hop to 3 out of 4 of its nearest 
neighbour sites only, dictated by an oriented network of constrained trajectories corresponding to the 
ensemble of classical Dirac strings [30] of connected dipoles [33] . The positively charged monopoles move 
in one sense along the network while the negative charges move in the opposite direction (see figure 4.b). 
The network is dynamically re-arranged through the evolution of the monopole configuration. In fact the 
characteristic time scale that we compare with experiment comes from the evolution of the network of 
Dirac strings rather than from the monopoles themselves. That is, the magnetic relaxation experiment 
of Snyder et. al. [45] ultimately measures spin dynamics, due monopole diffusion, rather than giving a 
direct measure of the monopoles themselves. We define the string network by an integer a = ±1 giving 
the orientation of the Dirac string along each bond of the diamond lattice. The Dirac string network is 
an image of the underlying spin structure for which we can define the auto-correlation function 



We start the simulation with the Dirac string network in an ordered configuration and let it evolve at 
T = 4 K until equilibrium is reached. A high density of monopoles is generated, which move following 
the criteria defined above and time step is defined as a sweep through N links of the monopole vacuum. 
We compute C(t) averaged over 200 simulations following the same procedure as for the calculation of 
spin-spin autocorrelation function in the DSL The Monte Carlo time was again scaled to the experimental 
time at T = 4 K and no other fitting parameters were used. 

5.3 Comparison with experiment 

In accordance with our analysis of figure 11, the characteristic time scale for relaxation of the Dirac 
string network does indeed reproduce the sharp spin freezing observed experimentally [27]. Results are 
shown in figure (12) in the interval below 4 K. Data from the Coulomb gas simulation is in very close 
agreement with the Monte Carlo data for the DSI, confirming that the low temperature dynamics of the 
DSI is driven by the dynamics of Coulombic quasi particles, so that a Coulomb gas in the grand canonical 
ensemble really is the correct description of the system at low temperature. Also shown in the figure is the 
modified Arrhenius formula 14 describing freezing in the NNSI. On this temperature scale the difference 
between data for the NNSI and that for the DSI or the Coulomb gas is considerable, confirming that 
the sharp, non-Arrhenius freezing at low temperature is due to the presence of long range interactions 
between quasi-particles. The most significant effect here is the gradual suppression of Debye screening 
below ~ 2 K, provoking the avalanche effect which accelerates the disappearance of the monopoles below 
this temperature. In addition we also expect the diffusive dynamics to be slowed by the formation of 
neutral pairs, bound by the Coulomb interaction. However, as the Coulomb interaction is not confining 
in three dimensions the monopole mobility should remain finite. Below 1.5 K the estimate of the time 
scale for the Coulomb gas lies slightly below that for the DSI. This small difference provides a measure 
of the approximation involved in passing from the spin system to the charged particles. 

We show simulation results using both our definitions of the chemical potential. The data sets are 
very close to each other, although a small systematic difference appears around 1 K, with data for the 
fixed chemical potential, /ii giving a slightly lower estimate for t. This small deviation seems consistent 
with the small evolution in n(T) vs T observed in figure (11). 

Finally, when looked at with this level of precision, all numerical data falls below the experimental 
data at low temperature. This small but systematic difference may indicate the limits of applicability 
of both the DSI and our model of stochastic dynamics although more experiments need to be done to 
quantify this remark. More recent analysis of the static properties of Dysprosium Titanate suggest that 
the DSI can be improved by including further neighbour cxchangc[56], which will change fine details of 
the dynamics. The approximation of having a single, temperature independent microscopic time scale, 
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To, will also have its limitations which will also be manifested in the fine detail of the difference between 
simulation and experiment. 




12 3 4 

Temperature T [K] 

Figure 12: Spin freezing in Dj^T^Oy: The experimental data (x) are from Snyder & al. [45]. The 
Arrhenius law (blue line) represents the free diffusion of topological defects in the nearest neighbour 
model. The relaxation time scale of the Dirac string network for the Coulomb gas has been obtained 
for fixed chemical potential (a) and with /i varying slowly to match the defect concentration in dipolar 
spin ice (■), whereas the dipolar spin ice relaxation time is given by •. As in figure 9 and as defined in 
the text, the bars do not represent the error bars on the value of r; rather they show the extreme values 
ti and T2 at each temperature. The system size is 8 x 12 3 = 13824 sites on the diamond lattice for the 
Coulomb gas and 8192 spins on the pyrochlore lattice for DSI simulations. 



6 Conclusion 

The recent discovery of monopole quasi-particles in spin ice models and materials provides a rare occasion 
to think outside the magnetic box, while studying a set of magnetic materials. The theoretical[30, 36], 
numerical[27, 54, 37] and experimental[17, 36, 41, 18] work so far presented in this field make this point 
very clear. Here, the collective quasi-particle excitations are localized in real space, rather than reciprocal 
space, despite the absence of quenched disorder. The quasi-particle vacuum has hidden internal degrees 
of freedom with U(l) gauge symmetry and thanks to the unique role played by the dipole interactions 
on the pyrochlore lattice structure, the excitations have both the topological and energetic properties of 
magnetic charges. The correct description at low temperature is therefore that of a Coulomb (lattice) 
gas in the grand canonical ensemble, whose properties are almost identical to those of its electronic 
equivalent, an electrolyte of positive and negatively charged ions whose dynamics can be considered as 
purely stochastic. 
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There are some exceptions to this analogy, notably spin ice materials cannot support a dc monopolc 
current in the steady state [27]. This is because of the dynamical constraint that positive and negative 
charges move in opposite directions over the underlying Dirac string network, changing the direction 
of the network as they pass. The system only posses a finite number of monopole trajectories in each 
direction before complete ordering of the network occurs. Or, in magnetic language, the sample can 
provide a finite number of spin flips in a given direction before becoming magnetically ordered. This fact 
seems difficult to reconcile with recent experiments [41], where the magnetic charge has been measured 
through application of Onsager's Wien effect [40] to Dysprosium Titanate. In an electrolyte, in the low 
charge concentration limit, the electric current density, a(E) increases linearly with electric field E, with 
the result that the ratio a(E)/a(0) is a universal number containing E, T, charge q and fundamental 
constants only. By treating muon relaxation rate in an equivalent way to the electric current, Bramwell 
et. al. have been able to measure the magnetic charge Q « 4.6^s A -1 in good agreement with the 
theoretical prediction of Castelnovo et. al. [30]. How can this be, given that no permanent current can 
exist ? One possible explanation is that, in weak fields and on short time scales, the constraints do not 
qualitatively affect the dynamics. One can define a drift length scale for field A, lx = ksT/QX, where 
X is the magnitude of a magnetic or electric field. Above this scale one should observe the deterministic 
drift of the ions in the electrolyte and the effects of the constraints in the magnetic problem. In the 
field and temperature range of the Wien effect experiment (T = 100 — 200 mK, B w 2 mT), Ib ~ 70 A. 
Given this large scale it seems reasonable that the magnetic and electric systems should show the same 
behaviour within the experimental set up. Much work is required here to understand these points in 
detail. 

Spin ice dynamics at low temperature should be compared with the dynamics of a spin glass. In 
the latter, the glassy behaviour is the result of a hierarchy of time scales, as manifested, for example in 
the phenomenological random energy model [57]. Here the evolution steps occur with an exponentially 
varying set of delay times. For spin ice, at least in zero field things are very different. We have modeled 
the dynamics of the NNSI using a single time scale associated with the quantum tunneling of an Ising 
spin ro- This scale dictates two processes; a hopping move of an existing deconfincd monopole by one 
lattice spacing, or the creation (annihilation) of a nearest neighbour pair of particles, whose characteristic 
time is dressed by the Boltzmann weight, r = r D exp(— (3AE). At low temperature this becomes the 
rate determining step, so that once created, the diffusing particles are actually extremely efficient at 
decorrelating the system. For example, from figure (7), the characteristic Arrhcnius time scale at 500 
mK is t ~ 1 sec, while the microscopic time scale, To is estimated to be about 2.5 miliseconds. That 
is, the system decorrelates in 400 microscopic time steps thanks to the free diffusion of the deconfincd 
particles. In the presence of Coulomb interactions, the decorrelation times are longer as the microscopic 
processes are dressed by a continuous range of energy scales coming from the Coulomb interaction, but 
the basic picture is provided by the simpler nearest neighbour model where free particles move on an 
essentially flat energy surface, making them efficient in decorrelating the system. 

In this paper we have interpreted a magnetic relaxation experiment in terms of monopole dynamics. 
A more fundamental approach is to design experiments that explicitly and directly test for monopole 
dynamics. Up to press the only experiments which attempt this are those using the Wien effect [41, 39], 
in which monopole mobility can be probed by exploiting the chemical equilibrium existing between bound 
and free particles in a field. An even more direct approach would be an equivalent of the "Stanford" 
superconducting coil experiment [58, 30], whose goal was to measure the passage of free Dirac monopoles 
through a superconducting coil. Such an experiment could in principle detect the passage of a single 
magnetic quasi-particle but, given that the charges have no mass and therefore have diffusive, rather 
than Newtonian dynamics, at least in the absence of an external field and that the monopole charge in 
Dysprosium Titanate is approximately 8000 times smaller than a Dirac Monopole [30], this experiment 
presents serious challenges. A starting point could be a many body equivalent. In fact, such an experiment 
has been set up in another context by Herrison and Ocio[59] to look at the violation of the fluctuation 
dissipation relation in spin glass systems. In this experiment the magnetic fluctuations are directly 
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measured and compared with magnetic response. The analysis presented in this paper suggests that the 
resulting magnetic noise in spin ice would be due to the diffusive motion of monopoles. Unfortunately, 
for the moment this set up is limited to the 4 K temperature range which remains in the high monopolc 
density regime. The goal for a monopole experiment would be first to descend below one Kelvin and then 
to reduce the scale of the sample sufficiently so that the resulting noise was due to a mesoscopic, or even 
microscopic number of quasi-particles. 

In conclusion, the discovery of magnetic monopoles in spin ice has exposed a rich vein of physics, 
in which a Coulomb gas, equivalent to a weak electrolyte appears at low temperature. This emergent 
'magnetolyte'[41] shows many of the static and dynamic properties of the electric system within a very 
idealized environment. Much experimental, theoretical and numerical work remains to be done in order 
to fully exploit it. 
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